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ABSTRACT 

For more than a decade now, it has been controversial whether or not the high rate of giant gravitational arcs and the largest ob- 
served Einstein radii are consistent with the standard cosmological model. Recent studies indicate that mergers provide an efficient 
mechanism to substantially increase the strong-lensing efficiency of individual clusters. Based on purely semi-analytic methods, we 
investigated the statistical impact of cluster mergers on the distribution of the largest Einstein radii and the optical depth for giant 
gravitational arcs of selected cluster samples. Analysing representative all-sky realizations of clusters at redshifts z < 1 and assuming 
a constant source redshift of z $ = 2.0, we find that mergers increase the number of Einstein radii above 10" (20") by ~ 35% (~ 55%). 
Exploiting the tight correlation between Einstein radii and lensing cross sections, we infer that the optical depth for giant gravitational 
arcs with a length-to- width ratio > 7.5 of those clusters with Einstein radii above 10" (20") increases by ~ 45% (~ 85%). Our findings 
suggest that cluster mergers significantly influence in particular the statistical lensing properties of the strongest gravitational lenses. 
We conclude that semi-analytic studies must inevitably take these events into account before questioning the standard cosmological 
model on the basis of the largest observed Einstein radii and the statistics of giant gravitational arcs. 
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1. Introduction 

Both the abundance of gravitational arcs and the distribution 
of Einst ein radii in galax y clusters are valuable c osmological 
probes (Bartelmann 2010 ). Therefore, the result of Bartelmann 
|et al.| ( |1998| ), who reported that we observe ten times as many 
giant gravitational arcs on the sky as theoretically expected, 
poses a serious challenge. Various aspects of strong gravitational 
lensing that could potentially mitigate the tension between the- 
ory and observations were studied in a lo ng series of subse- 
quent works (see e.g. |Fedeli et al.| ( [2008| ) or |B artelmann| ( |20 10| 
Sect. 5.2 ) for summaries, and references therein). In particular, 
|Torri et aL] d2004 ) analysed numerically simulated mergers of 



galaxy clusters. These authors found that mergers substantially 
change the shape of the critical curves and can boost a clus- 
ter's efficiency to p roduce giant arcs by an order of magnitude. 
|Fedeli et aklp006| ) employed semi-analytic methods to estimate 
that cluster mergers approximately double the statistical strong- 
lensing efficiency of clusters at redshifts z > 0.5. Fedeli and 
coworkers argued that mergers might possibly explain the ex- 
cess of gravitational arcs in obser ved galaxy cluster s at moder- 
ate and high redshifts. However, Fedeli et ak| p006 ) made sev- 



eral simplifying assumptions that we revise here: First, galaxy 
clusters were described by elliptically distorted spherical lens 
mod els instead of adopting more realist ic triaxial density pro- 
files dling & Suto 2002; Qgu ri et al.|20 03 ). This approximation 
reduces the required computing time substantially, since the cal- 
culation of deflection angles for triaxial density profiles involves 



numerical integrations ( |Schramm|1990b| ), while simple analytic 
expressions exist in the cas e of elliptically distorted density pro- 
files ([Schneider et al.|1992| ). Second, all mergers were simulated 
with a fixed direction of motion and relative orientation of the 
merging clu sters , neglecting two important degrees of freedom 
(see Section \5?2\ . 

Recent studies indicate that the distribution of Einstein radii 
might also be in conflict with theory. More precisely, the largest 
observed Einstein radii (e.g. |Halkola et al.||2008 |Umetsu &| 
[Broad hurst 2008; Zitri n et al.||20TlJ " were claimed to exceed 
the maximum possible expectations of the standard cosmolog- 
ical model ([B roadhurst & Barkana 2008; Oguri & Blandford 
2009; Meneghetti et al. 201 1|). These conclusions were drawn 
by either comparing the largest observed Einstein radii to those 
found in numerical simulations or by semi-analytically estimat- 
ing the probability of finding the strongest observed lens systems 
in a ACDM universe. While studies of Einstein radii in numer- 
ical simulations are probably most realistic, they always suffer 
from a limited sample size. The simulated boxes might simply be 
too small to contain a sufficient number of extraordinarily strong 
gravitational lenses, which forbids solid statistical conclusions. 
In a follow-up paper, we will show that this limitation is indeed 
decisive in the context of extreme value statistics ( Waizmann 
|et al.|2012| ). Semi-analytic methods - admittedly based on a set 
of simplifying assumptions - can overcome this limitation be- 
cause they are computationally less demanding and hence can 
be used to analyse large samples of particularly strong gravita- 
tional lenses within a comparably short time. However, we note 
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that so far all semi-analytic studies of cosmological distributions 
of Einstein radii have only considered samples of isolated galaxy 
clusters. One important goal of this paper is to extend these pre- 
vious approaches and to present a new semi-analytic method for 
studying distributions of Einstein radii that incorporates the im- 
pact of cluster mergers. 



New findings of Meneghetti et al. ( 201 1 ) suggest that the ex- 
cess of giant arcs and the too large Einstein radii are closely re- 
lated. Analysing selected samples of strong gravitational lense s 
in the MareNostrum simulation (Gottlober & Yepes 2007), 
Meneghetti and collaborators discovered a remarkably tight cor- 
relation between lensing cross sections and Einstein radii of 
cluster-sized dark matter haloes. This correlation plays an im- 
portant role for the structure and line of reasoning of this work. 

This paper is structured as follows. Section [2] summarizes 
an analytic model of triaxial gravitational lenses. In Section [3] 
we compare three different methods for computing lensing cross 
sections for giant gravitational arcs and show how the alternative 
algorithms need to be adjusted so that they yield equally reliable 
results. Having identified a reliable measure of lensing cross sec- 
tions, we first investigate the correlation between Einstein radii 
and lensing cross sections in Section [4] and then study how the 
correlation evolves during cluster mergers in Section [5] Finally, 
in Section [6] we introduce a semi-analytic method for comput- 
ing cosmological distributions of Einstein radii that properly in- 
cludes cluster mergers. Our conclusions are presented in Section 

m 

In this paper, we consider two different sets of cosmo- 
logical parameters obtained from the Wilkinson Microwave 
Anisotropy Probe (WMAP). These are the b est-fitting pa- 
rameters from the WMAP one-year (WMAP1; Spergel et al. 
12003) , (QmcQAo.nbo, A, cr 8 ) = (0.7, 3,0.045,0.7,0 . 9) and 
from the seven-year data (WMAP7; Larson et al. |2011 
(Q m o,n AO ,Q b o,/*,0- 8 ) = (0.727,0.273,0.0455,0.704,0.811). 



2. Triaxial gravitational lenses 

2.1. The density profile of triaxial dark matter haloes 



Jing & Suto ( 2002 , hereafter JS02) performed a combined anal- 
ysis of triaxiality in twelve high-resolution dark matter halo 
simulations to study the detailed shape of single haloes as 
well as haloes embedded in large cosmological simulations, to 
gather statistical information about halo profiles. These authors 
found that the universal, spherical density profile discovered by 
Nav arro et al.| ( |1996| hereafter NFW) can be generalized to a 
triaxial model and showed that this generalization significantly 
improves the fit to simulated haloes. Moreover, by analysing 
large cluster populations in their cosmological simulations, JS02 
derived probability density functions for the profile concentra- 
tions and axis ratios of triaxial dark matter haloes. The statis- 
tical description provided in JS02 allows constructing random 
catalogues of triaxial dark matter haloes that resemble realistic 
cosmological populations in numerical simulations. 

|Jing & Suto] ( |2002| ) parametrized the spatial density profile 
of a triaxial dark matter halo by means of Cartesian coordinates 
x' = (x\y',z') in the principal axis frame. Using this parametri- 
sation, they proposed a generalization of the NFW density pro- 
file 
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Here, z denotes the halo's redshift, Rq is the scale radius (cf. Eq. 
fT3])), S CQ is the characteristic density (cf. Eq. ([14])) and p C nt(z) 
denotes the critical density of the universe. The exact numeri- 
cal value of the inner slope a of the density profile is still be- 
ing discussed. While NFW originally proposed a - 1.0, other 
authors argued that steeper profiles with values up to a = 1.5 
provide a better fit to observations and numerical simulations 
(|Moore et al.|1999[|Jing & Suto|2000||Power et al. 12003} |Navarro| 
et al.|2004||Limousin et al.|2008| ). In contrast, recent combined 
strong- and weak-lensing analyses of selected cluster s indicate 
shallower mass profiles with inner slopes a < 1 (see Newman 



|et al.|20lT| for instance). Following Qgu ri et aL| ( |2003| ), we con- 
sider both a = 1 and a = 1.5 to cover a broad range of the 
predicted values and discuss some consequences of varying in- 
ner slopes on strong-lensing statistics. 

To draw a random triaxial dark matter halo of virial mass M 
at redshift z, we first sample its axis ratios using the empirically 
derived probability density functions from JS02, 
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Here, M* is the characteristic nonlinear mass scale where 
cr(M*,z) = S c (z). According to JS02, the best-fitting parameter 
for the width of the axis ratio distribution p(a/c) is cr s = 0.113. 
The concentration c e is defined by 



C e = 



Rz 



(7) 



where R Q is determined such that the mean density within the el- 
lipsoid of the major axis radius R e equals A e (z)£Xz)p C ritfe), where 



A e (z) = 5A vir (z) 



ab 



0.75 
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A v i r (z) is the overdensity of obj ects virialized at redshift z, which 
we approximate according to Nakam ura & Suto| ([1997 ). JS02 
found a log-normal distribution for the concentration, 



P (c Q ) = 
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- exp < - 
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with a dispersion of <x c = 0.3. Following Qg uri et aL] p003| ), 
we include a correlation between the axis ratio a/c and the mean 
concentration, 



c q — fc 



A vir (z c ) / l +z, 
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In Eq. fTTT), we adopted a correction introduced by Qg uri & 



Blandford (2009), forcing f c > 0.3 to avoid unrealistically low 
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concentrations for particularly low axis ratios a sc . Additionally, 
as noted earlier by e.g. |Oguri & Keeton (2004), triaxial haloes 
with particularly low axis ratios a sc (and hence also low concen- 
trations c Q ) are highly elongated objects whose lens potential is 
dominated by masses well outside the virial radius. While pre- 
vious studies (e.g. |Baltz et al.|[2009| |Oguri & Blandford||2509] ) 
tried to avoid these unrealistic scenarios by simply truncating 
the density profile (cf. Eq. ([l}) outside the virial radius, we alter- 
natively force all sampled axis ratios a sc to lie within the range of 
two standard deviations. Since this manipulation may affect the 
statistics of the largest predicted Einstein radii, we will investi- 
gate^ the consequences in our follow-up work ( Waizm ann et al. 
|2012| ). The numerical value of the free parameter A e depends on 
the underlying cosmology. Throughout, we set A e = 1.1, which 
was proposed by JS02 for a standard ACDM model. The above 
expressions only apply to an inner slope of a = 1 .0. In the case of 
a - 1.5, we used the simple relation c Q (a = 1.5) = 0.5 x c Q (a - 
1.0) ( [Keeton & Madau|p00l) ; JS02). z c denotes the typical col- 
lapse redshift of a dark matter halo of mass M and is computed 
using the complementary error function, 



erfc < 



{ V2 [cr 2 (fM) - cr 2 (M)] J 2 



(12) 



Here, <x 2 (M) is the top-hat smoothed variance of the power- 
spectrum extrapolated to the redshift z = and oj(z) = 
S c (z)/D(z), where D(z) is the linear growth factor. Following 
JS02, we adopted / = 0.01. The typical collapse redshift can 
be derived within the framework of extended Press-Schechter 
theory ( [Lacey & Cole|1993] ). It corresponds to the typical time 
when the most massive progenitor of a dark matter halo con- 
tained the mass fraction fM. However, note that JS02 defined 
the typical collapse redshift such that it does not depend on the 
halo's actual redshift. 

All other profile parameters can be inferred from the axis 
ratios and the concentration. Using an empirical relation (found 
by JS02) between R e and the spherical virial radius r v i r , R e /r Y i Y « 
0.45, the scale radius Rq is given by 



R = 0.45 
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Finally, the characteristic density S CQ of the density profile ([]} is 
defined in terms of the concentration, 
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2.2. Basics of gravitational leasing 

Since the theory of gravitational lensing is well known, we ab- 
stain from giving a thorough introduction but instead refer to 
( |2010| ) for a recent review. We only summarize the 



Bartelmann 



quantities that are particularly important here. 

Working in the thin- screen approximation, gravitational 
lensing essentially reduces to a mapping between the lens plane 



and the source plane. The lens mapping is given by the lens equa- 
tion, 

/3 = 0-a(0) , (16) 

which describes the relation between the position ft of a source, 
the position of the lensed image and the reduced deflection an- 
gle a (0). The magnification introduced by gravitational lensing 
is given by 

(i7> 

where A t and A r denote the tangential and the radial eigenvalue 
of the lens mapping. Points on the lens plane where at least one 
eigenvalue vanishes are called critical points. The set of critical 
points forms closed curves called critical curves. Their images 
on the source plane are called caustics. Equation ( [T7] ) indicates 
that sources located close to caustics are strongly magnified by 
the lens mapping. Their images - such as giant arcs - appear in 
the regions close to the critical curves in the lens plane. 



2.3. Gravitational lensing by triaxial dark matter haloes 

Equation ([T]) defines the density profile of triaxial dark matter 
haloes in their principal axis frame. However, their lensing prop- 
erties are determined by their projected surface mass density in 
the observer's frame. Since observed haloes may have arbitrary 
orientations and the required line-of-sight projection requires us 
to express the spatial density profile in terms of coordinates with 
respect to the observer, we first introduce an appropriate coordi- 
nate transformation. 

Adopting the notation of Eq. ([T]), x' = (x',y',z') denote the 
Cartesian coordinates in the principal axis frame of the halo, 
so that the z'-axis lies along the major axis of the ellipsoid. 
Moreover, we assume that the z-axis is aligned with the line- 
of-sight direction of the observer. Then, we can locally construct 
another Cartesian frame for the observer's coordinate system, 
denoted by x = (x,y,z). The origins of both coordinate sys- 
tems are placed at the centre of the halo. |Oguri etaL ( |2003 ) 
exclusively investigated the lensing properties of individual tri- 
axial haloes and hence could fix the relative orientation of the 
projected surface mass density profile in the x-y plane at will. 
This freedom allowed them to use a simplified rotation matrix to 
parametriz e the coordinate transformation between both frames. 
In Section 5.2 we show that the relative orientation of the pro- 
jected mass density profiles of merging clusters plays an impor- 
tant role and accordingly we cannot adopt this practical sim- 
plification. Instead we have to introduce a general coordinate 
transformation between both frames. To this end, we introduce 
the three Euler angles (ifs,6,(p) and adopt the so-called z-x'-z!'- 
convention illustrated in Fig. [I] The halo's coordinate system is 
rotated about the z-axis by the angle 0, then about its new x'-axis 
by the angle 6, and finally about the new z"-axis by the angle 0. 
The coordinate transformation is given by 



x' = Rx , 



(18) 



where the rotation matrix R is the product of the three rotation 
matrices about the axes. Using this coordinate transformation, 
we can express the density profiles of arbitrarily oriented triaxial 
dark matter haloes in terms of the observer's coordinates. From 
here on, the derivation of all lensing pro perties of triaxial dark 
matter haloes follows Qgu ri et al. ( 2003| ), hence we refer the 
reader to their work. However, since we had to introduce a more 
general coordinate transformation, some algebraic coefficients 
are more complicated in our approach. We provide the necessary 
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perimeter found from boundary points 
• image points 



length correction 



Fig. I: Visualisation of the Euler angles in the z-x'-t!' conven- 
tion, which are used to parametrize the coordinate transforma- 
tion between the observer's frame and the principal axis frame 
of the triaxial dark matter halo (figure adapted from Wikipedia 
(Lionel Brits)). 



expressions in Appendix [A] Moreover, we present a practical 
method to speed up the computation of deflection angles by a 
factor of ~ 10 in Appendix [B] 

In the following sections, we analyse strong-lensing proper- 
ties of randomly drawn cosmological populations of triaxial dark 
matter haloes. Therefore we also need to sample random halo 
orientations. For that purpose, special care must be taken, since 
we cannot simply uniformly distribute the three angles (if/, 0, (p) 
in the corresponding an gle bins. In stead, we have to apply the 
correct Haar measure (Haar 1933| ). Only if/ and can be uni- 
formly distributed in the range [0, 2n]. To sample 6, we first draw 
a uniformly distributed random number r in the range [0, 1) and 
then compute 

6 = arccos(l - 2r) . (19) 



3. Computation of lensing cross sections 

3.1. Definition of the lensing cross section 

The strong-lensing efficiency of an arbitrary mass distribution 
is typically measured by the so-called lensing cross section cr p , 
which is defined as the area of the region on the source plane 
where a source with given characteristics (morphology, orien- 
tation, surface brightness profile, etc.) has to lie in order to pro- 
duce at least one image with properties p. This definition is quite 
general and needs to be additionally specified for the following 
applications. Throughout, we exclusively computed the lensing 
cross section cr 7 5 , which measures the efficiency of a gravita- 
tional lens to produce highly elongated images (i.e. giant arcs) 
with length-to- width ratios exceeding 7.5. For determining this 
area, each of the three following methods places random sources 
in the relevant regions of the source plane and analyses their 
le nsed images. As was shown by e.g. Bartelma nn et al.| ( [l995| ) 



Oguri (2002), the size and the morphology of the background 



or 

sources influence the resulting lensing cross section for gravi- 
tational arcs significantly, therefore a precise definition of these 
characteristics must be given. Following the common convention 
of previous works (e.g. Bartelmann et al. 1995 ; Meneghetti et al.| 
[20001 [2003} |Puchwein et al.||2005| ), we model the background 
sources as elliptical galaxies with constant surface brightness. 




Fig. 2: Illustration of the image fitting algorithm and the length- 
and perimeter-correction of the ray-tracing method A (figure 
from |Puchwein et aL ( 2005| )). 



Their size is chosen such that it equals the size of a circular 
source with a radius of r = 0.5". Their orientation (in the source 
plane) and minor- to-major axis ratios are randomly drawn from 
the ranges [0,7r] and [0.5, 1], respectively. 

3.2. Ray-tracing method A 

The first method for computing lensing cross sections was in- 
troduced by Miralda-Escude (1993), adopted to nonanalytic lens 
models by Bartelmann & Weiss ( 19 94]) and further extended by 
subsequent works (Menegh etti et al.| 2000, 2003; Puch wem et al.| 
2005). Elaborate descriptions of this method including all tech- 
nical details are given there, hence we restrict ourselves to a brief 
summary of the aspects relevant for the following discussion. 

First, the lens plane is covered with a Cartesian grid. At each 
grid point, the light deflection from gravitational lensing is com- 
puted to end up with a map of deflection angles. The resolution 
of this map must be sufficiently high so that the simulated images 
are resolved accurately and pixelisation effects are negligible 
(see below). Second, the source plane is covered with a Cartesian 
grid whose resolution is adaptively refined in regions close to 
the caustics. Next, the algorithm loops over each point of the re- 
fined grid and repeats the following set of operations. A random 
elliptical source is placed on the current grid point. Then, the 
images of this source are simulated using ray-tracing. A bundle 
of light rays starting at the observer is traced through the lens 
plane and mapped (by means of the lens equation) to the source 
plane. If a ray hits the source, the corresponding pixel on the lens 
plane is marked as an image point. After that procedure, the dis- 
tinct resulting images are identified using a standard connected- 
component labelling algorithm and are separately analysed by 
a simple image-fitting algorithm, which is illustrated in Fig. [5] 
This algorithm sequentially determines three characteristic im- 
age points: 

1. Point (a), which represents the arc centre and is given by the 
point whose image on the source plane falls nearest to the 
source centre. 

2. Point (b) on the boundary, which is farthest away from the 
arc centre (a). 

3. Point (c) on the boundary, which is farthest away from (b). 

After that, a circle is fitted through the points (a), (b), (c) 
and the length of the arc connecting these three points is pre- 
liminarily set as the image length L. Next, the perimeter is mea- 
sured by walking along the ordered boundary points and sum- 
ming up their mutual distances. Since the images are simulated 
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using d iscrete grids, pixelisati on effects have to be taken into ac- 
count. Puchwein et al. ( 2QQ5| ) correctly noted that cells on the 
lens plane are only classified as belonging to the image if their 
centres fall within the source. Thus, boundary points are on aver- 
age half a grid constant farther inside the image. To compensate 
for this effect, four grid constants must be added to the mea- 
sured perimeter and the measured length must be elongated by 
one grid constant. The area A of the image is found by simply 
counting the number of image pixels (blue points in Figure [2} 
and summing up their corresponding cell sizes on the lens plane. 
Next, four different figures are fitted to the image: 

1. A rectangle with area A, length L, and setting the width to 
W = AIL. 

2. An ellipse with area A and assuming that the length of the 
major axis equals the measured image length L. The image 
width W is then defined as the length of the minor axis, W = 

A/(7TL). 

3. A ring with outer radius r out = L/2 and choosing the inner 
radius r m such that the area of the ring equals the measured 
area A of the image. The image length is then redefined as 
L = n(r out + ri n ) and the width is set to the width of the ring, 
W = r out - n n . 

4. A circle with radius r = ^A/n and hence defining r = L = 
W. 

After this fitting routine, the algorithm determines the best 
approximation of the real image shape by singling out the fig- 
ure whose circumference agrees best with the measured image 
perimeter. If the length-to-width ratio of the chosen figure is 
higher than 7.5, the corresponding cell size on the source plane 
is added to the lensing cross section cr 7 5 . If the length-to- width 
ratios of multiple images of the considered source exceed that 
threshold, the cell size is multiplied by the corresponding image 
number. 

Performing the above procedure for all points of the adap- 
tive grid on the source plane yields an accurate estimate of the 
lensing cross section cr 7 5 . 

3.3. Ray -tracing method B 

During a series of tests with method A, we observed th at remov- 
ing the lengt h- and perimeter-correction introduced by Puchw ein 
et al.| ( |2005| ) on average decreases the computed lensing cross 
sections by roughly 25%, indicating that (unavoidable) errors in 
the length- and perimeter-measurements also have a significant 
impact on the results. Moreover, realistic shapes of gravitational 
arcs are rather irregular and only barely agree with the fitted, ide- 
alized geometric figures. To cross-check the results of method 
A, we developed a new ray-tracing method that employs an al- 
ternative algorithm to determine length-to-width ratios of lensed 
images. As in method A, randomly drawn elliptical sources are 
placed on an adaptively refined grid on the source plane. Again, 
the corresponding image configurations are found by tracing 
bundles of light rays through the lens plane. To identify the dis- 
tinct lensed images, we employ the component-labelling algo- 
rithm proposed by Chan g et aL] ([2004 ), since it returns a clock- 
wise ordered list of boundary points by construction. This prop- 
erty is important for our alternative length-to-width measure- 
ment, which is illustrated in Fig. [3] As with method A, we first 
identify the characteristic image points (a), (b), and (c) and fit a 
circle through them. Next, we determine the length of the (green) 
arc and correct for pixelisation effects by adding one grid con- 
stant. After that, we use points (b) and (c) to split the boundary 



length 




Fig. 3: Illustration of the new algorithm to measure length-to- 
width ratios of gravitational arcs. The image length is measured 
by fitting a circular arc (green arc) through the points (a), (b) 
and (c). Instead of fitting idealized, geometric figures to infer 
the width of the image, the average width is measured along the 
arc. 



into two parts. The first part consists of pixels connecting (b) and 
(c) in clockwise direction (upper part of the boundary) and the 
second group of pixels connects both points in anti-clockwise 
direction (lower part of the boundary). This separation of bound- 
ary points allows us to measure the average width of the image. 
For this, we first subdivide the angle that is enclosed by the arc 
connecting (b) and (c) into small bins. Then, we consider each 
angular bin separately, scan both boundaries for pixels that lie 
in the current range and use them to determine the image width 
at the considered position. Figure [3] visualizes this procedure for 
two specific angular bins along the arc. Finally, we define the av- 
erage of all measured widths as the width of the image. Again, 
we correct for pixelisation effects by adding one grid constant. 
In the same way as in method A, we determine the number of 
images whose length-to- width ratio is higher than 7.5 for each 
source. We multiply this number with the corresponding cell size 
on the source plane and add this area to the final lensing cross 
section cr 7 5 . 

When we had completed this new algorithm, we realized that 
Menegh etti et al.| ( [2008| ) used the same approach to measure the 
length-to- width ratio of arcs in realistic observations. But instead 
of averaging the widths of the distinct angle bins, they computed 
the median in order to be less sensitive to observational noise, 
which produces extreme outliers. Because we did not experi- 
ence these problems in our purely theoretical computations, the 
simple average is expected to yield similar results. 

3.4. Semi-analytic method 

These ray-tracing simulations have a significant drawback: they 
can become computationally quite expensive. As already men- 
tioned, the ray-tracing part requires computing a map of deflec- 
tion angles. The resolution of this map must be sufficiently high 
so that pixelisation effect^] do not significantly falsify the mea- 
sured length-to- width ratios. Performing a series of convergence 
tests, we found that the images of unmagnified circular sources 
with radius r = 0.5" should at least be represented by ~ 250 
pixels on the lens plane, which amounts to a grid constant of 
~ 0.056". The deflection angle map needs to cover the entire 
region of the lens plane in which gravitational arcs can occur. 



1 The perimeter measurement is particularly sensitive to pixelisation 
effects. When images have an unfortunate orientation with respect to the 
Cartesian grid on the lens plane, their boundary, which should normally 
be smooth, is fairly rugged owing to their discrete representation on the 
grid. Thus, by summing the mutual distances between the ordered list 
of boundary points, the measured perimeter overestimates the correct 
circumference. 
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Hence, for particularly strong gravitational lenses, the required 
side lengths can easily reach values of ~ 200" or more, which 
in turn lead to Cartesian grids with roughly 3500 x 3500 pix- 
els or more. For triaxial gravitational lenses, the calculation of 
deflection angles involves numerical integrations and hence the 
computation of the deflection angle map can become rather ex- 
pensive. Moreover, we note that the component-labelling algo- 
rithms, which are the most expensive operation of the ray-tracing 
methods, require more and more time as the resolution of the de- 
flection angle map is increased because images are represented 
by ever larger ensembles of pixels. Finally, the resolution of the 
adaptive grid in the source plane must also be sufficiently high 
for a precise estimate of the lensing cross section. Consequently, 
many random sources must be placed and their images need to 
be simulated. 

If arc statistics is meant to be used as a cosmological 
probe, lensing cross sections of large cluster samples need to 
be computed. In addition, these computations need to be re- 
peated for different sets of cosmological parameters. Given the 
above requirements, this task quickly becomes computationally 
extremely demanding. To overcome this problem, |Fedeli et"aL] 
(2006) developed a fast and elegant method to calculate lensing 
cross sections semi- analytically. They found that their alterna- 
tive method is faster than ray-tracing simulations by a factor of 
~ 30 and yields equally reliable results. Since we aim to compare 
their method to the previously discussed ray-tracing simulations, 
a basic understanding of their new approach is required. Again, 
we only briefly outline the principal ideas and the most impor- 
tant formulae. We refer the reader to jFedeli et"aL] ( |2006| ) for a 
thorough introduction including all technical details. 

We note in the beginning that the intrinsic ellipticity of 
sources has a significant imp act on the lensing cross section 
for gian t gravitational arcs (cf. Bartelma nn et al.|1995| ). Keeton 



( 2001b) developed a simple and elegant formalism to treat this 



aspect mathematically, which we adopt as described by Fedeli 



|et al.| ([2006 ). Without going into detail here, we assume that 
we possess an appropriate method for computing the average 
of the eigenvalue ratio (q\) (q\ = max[\A t / A Y \,\A Y /A t \]) over ar- 
bitrary images at arbitrary positions in the lens plane. Then, the 
area of the stripe in the lens plane where gravitational arcs with a 
length-to- width ratio equal to or higher than 7.5 occur is simply 
given by 

B t = J H[(q ohs )(0) -7.5] d 2 6 , (20) 

where H (x) is the Heaviside step function and the integration is 
carried out over the lens plane. The notation (g bs) indicates that 
we use the average eigenvalue ratio (q\) for computing the ob- 
served axis ratio. Using the lens equation, B\ can be mapped to 
an area B s on the source plane, which is by definition the lensing 
cross section cr 7 5 we are searching for. Taking the magnifica- 
tion caused by the lens mapping into account, the lensing cross 
section is thus given by 

(Ti.5 = J H [<?obs>(0) " 7.5] | det A (0) \ d 2 6 , (21) 

where detA(0) denotes the Jacobian determinant of the lens 
mapping. Effectively, this elegant approximation reduces the 
lensing cross section computation to a simple area integral, 
which can be evaluated numerically. 



3.5. Comparison of the different methods 

We are now in the position to carry out the comparison between 
the three different methods to check their reliability. To this end, 
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Fig. 4: Lensing cross sections cr 7 5 in the range of 0-400 arcsec 2 
according to method A (blue crosses), the semi-analytic method 
(turquoise crosses) and method A fitting ellipses only (yellow 
triangles) as a function of the lensing cross sections cr 7 5 com- 
puted with method B for a sample of 300 massive dark matter 
haloes (see text). The red dashed line indicates the bisector of 
the cr 7 5 - cr 7 5 plane. The violet dashed line represents the best 
linear fit of the semi-analytic lensing cross sections. 



we sampled 300 random triaxial haloes placed at redshift z\ - 
0.5 and computed their lensing cross sections 0-7.5 using each 
of the three algorithms. The halo masses were logarithmically 
uniformly distributed within (5 x 10 14 - 2 x 10 15 ) M /z _1 , and we 
adopted a - 1.0 for the inner slope of their density profile. The 
source plane was placed at redshift z s = 2.0. 

Figure[4]shows the lensing cross sections 0-7.5 computed with 
ray-tracing method A and the semi-analytic method as a func- 
tion of the lensing cross sections 0-7.5 computed with ray-tracing 
method B. For reasons to be clarified below, we also plot the re- 
sults according to method A but fitting only ellipses to the simu- 
lated images. 

Obviously, both ray-tracing methods agree excellently in the 
range 0-7.5 < 225 arcsec 2 , although they employ inherently dif- 
ferent algorithms to measure length-to- width ratios of lensed im- 
ages. For lensing cross sections 0-7.5 > 225 arcsec 2 , the results 
computed with method A tend to be slightly higher. However, 
we solely attribute these deviations to technical limitations. The 
implementation of method A was written in Fortran-77 and we 
experienced problems when allocating large amounts of mem- 
ory. Accordingly, we had to limit the maximum possible grid di- 
mension for the deflection angle maps to 2048 x 2048 pixels. For 
the largest shown lensing cross sections (0-7.5 > 225 arcsec 2 ), the 
considered field sizes might simply be too large for the highest 
possible resolution, leading to coarsely resolved images. Hence 
the errors due to pixelisation effects increase and account for the 
differences between both ray-tracing simulations. Given these 
explanations, we conclude that both ray-tracing methods yield 
equally reliable results. 

In contrast to Fe deli et aL| ( |2006| ), however, we cannot re- 
produce the good agreement between the semi-analytic method 
and the ray-tracing simulations. Instead, the semi-analytic lens- 
ing cross sections are systematically too small. In order to quan- 
tify the deviations, we computed the best linear fit: 

o- 7 .5,sa = (0-583 ± 0.002)0-7.53 - (0.931 ± 0.287) . (22) 
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Here, o"7.5 ? sa and o~ 7.53 denote the lensing cross sections com- 
puted with the semi-analytic method and method B, respectively. 
In Fig. |4| the above linear relation ( [22] ) is indicated by the vio- 
let dashed line. Evidently the scatter around the linear relation 
is small, indicating a strong correlation between both quantities. 
The small fluctuations with respect to the linear fit can be ex- 
plained by different realisations of the randomly drawn source 
populations: as shown by Keeton (2001b), the observed length- 
to- width ratios of lensed images depend on the sources' intrinsic 
ellipticities and relative orientations with respect to the eigen- 
vectors of the magnification tensor. Since these properties are 
randomly drawn when placing the sources on the adaptive grid, 
some cells only contribute to the lensing cross section if the 
source placed on top of them has favourable properties. Thus 
the lensing cross sections experience statistical fluctuations. 

The above observations suggest that there must be a fun- 
damental difference between the semi-analytic method and the 
ray-tracing simulations. As already mentioned, we additionally 
computed the lensing cross sections with a slightly modified ver- 
sion of method A. Instead of probing four different geometri- 
cal shapes, we fitted the simulated images with ellipses only. 
Obviously, this modified ray-tracing algorithm yields results that 
agree well with the semi-analytic method. This information turns 
out to be useful in the next section. 




50 100 150 200 250 300 350 400 
(77.5 [arcsec 2 ] (method B) 

Fig. 5: Lensing cross sections 0-7.5 in the range of 0-400 arcsec 2 
computed with the modified semi-analytic method as a func- 
tion of the lensing cross sections 0-7.5 calculated with ray-tracing 
method B. The red dashed line indicates the bisector of the 0-7.5- 
0-7.5 plane. The lensing cross sections were computed for a sam- 
ple of 300 massive triaxial dark matter haloes (cf. Section [33]). 



3.6. Reconciling the semi-analytic method with the 
ray-tracing simulations 

The findings of the previous sections suggest a simple way for 
reconciling the semi-analytic method with the ray-tracing sim- 
ulations. We recall that we found excellent agreement between 
both ray-tracing methods although their algorithms for determin- 
ing the width of lensed images differ substantially. However, fur- 
ther analysis reveals that method A fits ~ 99% of the detected 
arcs with rectangles, which effectively equals method B's mea- 
suring the average image widths without imposing additional 
geometrical constraints. Method A yields almost identical re- 
sults to the semi-analytic method if only ellipses are fitted to 
the simulated images. These findings simply confirm that the 
semi-analytic method works as it should: By assuming that the 
averaged observable axis ratios (g bs) equal the length-to- width 
ratios of images, the method implicitly assumes that the lensed 
images roughly look like ellipses. Hence the good agreement 
with the modified version of method A is rather reassuring than 
indicating a failure of the semi-analytic approximations. 

We thus find a simple solution reconciling the semi- analytic 
method with the ray-tracing simulations. We therefore consider 
an arbitrary image with area A and length L. Fitting a rectangle, 
we first obtain the width W = A/L, and thus the length-to- width 
ratio is given by 

L L 2 

— = — • (23) 
W A v } 

Alternatively, we can also fit an ellipse. Using the following re- 
lation between the area of an ellipse, its length and its width, 



L W 
A = 71— — 

2 2 

we obtain the length-to- width ratio 



L 
W 



ttL 2 
4 A 



(24) 



(25) 



Obviously, if we fix the area A and the length L (as both ray- 
tracing methods do), rectangles have higher length-to- width ra- 



tios than ellipses and the difference is simply given by a fac- 
tor of 71 1 A. That is precisely the reason why lensing cross sec- 
tions decrease by ~ 40% if only ellipses are fitted to the sim- 
ulated images. The measured length-to- width ratios are lower 
and therefore fewer images are recognized as giant arcs. On the 
other hand, this suggests an easy solution for adapting the semi- 
analytic method to the ray-tracing simulations. Recalling that 
the majority of arcs seems to be better approximated by rectan- 
gles, we can simply multiply the calculated observed axis ratio 
#obs by A 1 71 to estimate the correct length-to- width ratio of most 
images. This simple multiplication is the mathematical equiva- 
lent of fitting images with rectangles instead of ellipses. To ver- 
ify this idea, we sampled another random population of triaxial 
haloes following the same recipe as in the previous section and 
computed their lensing cross sections with the modified semi- 
analytic method and method B. Figure [5] clearly indicates that 
the agreement between both methods is now reasonably good 
and the deviations typically remain well below 15%. 

We therefore note that the significant and systematic devia- 
tions between the semi-analytic method and the ray-tracing sim- 
ulations are simply caused by different definitions of length-to- 
width ratios of images. Gravitational arcs have irregular shapes 
and barely agree with idealized, geometrical figures, leaving no 
natural definition of a length-to- width ratio. Since Fig. [4] and 
Eqs. ([23])-([25]) clearly show that the actual definition has a signif- 
icant impact on the final lensing cross section, it is particularly 
important to state precisely which convention is used to yield 
comparable results. Moreover, the semi-analytic method can be 
reconciled with the ray-tracing simulations by introducing a sim- 
ple correction, so that it still provides a fast, reliable alternative 
for computationally demanding problems. Nevertheless, for the 
remainder of this paper, we always use ray-tracing method B to 
compute of lensing cross sections. 
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4. The correlation between Einstein radii and 
lensing cross sections 

4.1. Two definitions of Einstein radii 

Einstein radii measure the size of the tangential critical curve. 
While a natural choice exists for the Einstein radius in the case 
of axially symmetric lenses, it is not obvious how this concept 
can be transferred to arbitrary lenses with irregular tangential 
critical curves. For this r eason various alternatives were pro- 
posed by several authors (Broadhurst & Barkana 2008; Oguri 
|& Blandford1[2009l |Zitrin et al.||2011[ |Meneghetti et al.||2011| ). 
The following two definitions turned out to be most useful for 
our purposes. 

The first definition is of statistical nature. Let t denote 
the set of tangential critical points. The median Einstein radius 
#E,med is defined as the median distance of these points with re- 
spect to the lens' centre C , 



#E,med = median j (6 ifX - 6 c , x f + [o Uy - O^f | 9 t e G t 



(26) 



This definition was introduced by [Meneghetti et al.| (2011 



who analysed strong-lensing properties of numerically simu- 
lated clusters with irregular shapes that also showed substruc- 
ture. The centres of those clusters were defined as the location 
of the maximum of the smoothed projected density distribution. 
We note that the median was chosen in order to be less sensitive 
to extreme outliers of irregular critical curves. 

The second definition is of geometrical nature. Let A denote 
the area enclosed by the tangential critical curve. Then, the ef- 
fective Einstein radius #E,eff is defined by 



?E,eff 



(27) 



such that a circle with radius #E, e ff has the area A. Evidently, 
the median Einstein radius and the effective Einstein radius both 
agree with the original definition for axially symmetric lenses. 

4.2. Correlation for numerical clusters 



Menegh etti et al. ( 2011| ) extracted a sample of massive clus- 
ters [M > 5 x 10 14 M o /z _1 ) at redshifts z\ > 0.5 from the 
MareNostrum Universe (Gottlo ber & Ye pes 2007) and per- 
formed a statistical analysis of their strong-lensing properties. 
Figure [6] shows parts of their results, indicating a tight corre- 
lation between the lensing cross sections cr 7 5 and the median 
Einstein radii #E,med of their cluster sample. Performing a lin- 
ear fit to all data points in the log(#E,med) _ log(cr 7 5) plane, they 
obtained 

log(cr 7 . 5 ) = (1.79 ± 0.04) log(% med ) - (5.16 ± 0.05) , (28) 

with a Pearson correlation coefficient of r = 0.94, confirming 
the tight correlation between both quantities. 

These findings are important for several reasons. First, as al- 
ready noted by Menegh e'tti et al.| ( |2011| ), they clearly indicate 
the strong connection between the problem of too large Einstein 
radii and the arc statistics problem. A particularly strong gravi- 
tational lens, whose Einstein radius is too large for the ACDM 
model, will most likely also have a lensing cross section that ex- 
ceeds the maximum theoretical expectations. Hence, if we were 
to observe too large Einstein radii, we would also expect an ex- 
cess of gravitational arcs. On the other hand, we note that the 




100 



E ,med [arcsec] 



Fig. 6: Strong correlation between the median Einstein radii 
#E,med and the strong-lensing cross sections cr 7 5 of a sample of 
massive clusters with masses M > 5 x 10 14 M Q h~ l and redshifts 
Z\ > 0.5 extracted from the MareNostrum Universe. The red 
solid line indicates the best linear fit relation between log(<x 7 .5) 
and log(#E,med) (figure from Meneghetti et al . I d20 1 1 b) . 



converse is not true and care must be taken, since both problems 
are not identical. Even if there were a way to explain the partic- 
ularly large Einstein radii within the ACDM model, this would 
not necessarily solve the arc statistics problem. While the dis- 
tribution of largest Einstein radii only tests the extreme cases, 
an arc statistic additionally measures the cumulative lensing ef- 
ficiency of the entire halo population, so that it is also sensi- 
tive to the pure abundance of relatively unspectacular lenses, for 
instance. Thus, if the observed universe simply contained more 
strong gravitational lenses than theoretically predicted, we could 
hypothetically still observe an excess of gravitational arcs with- 
out having the problem of too large Einstein radii. 

Secondly, the correlation discovered might enable us to com- 
pute arc statistics for certain halo samples following a new ap- 
proach. Even by means of the fast semi-analytic method, the cal- 
culation of lensing cross sections for large halo samples is a com- 
putationally demanding problem, since complete maps of deflec- 
tion angles have to be computed for each individual lens. In com- 
parison, the calculation of Einstein radii can be fast because only 
few deflecti on a ngles near critical curves need to be computed 
(see Section [6^2]). Hence, the idea is to sample cosmological dis- 
tributions of dark matter haloes, compute their Einstein radii, and 
finally infer their lensing cross sections by means of the correla- 
tion. To thi s end, we need to investig ate if we can reproduce the 
findings of Meneghe tti et al.| ( |2011| ) with analytic mass profiles 
for triaxial haloes. 



4.3. Correlation for clusters with analytic mass profile 



We mimicked the cluster sample from Meneghetti et al. 



( 201 1) by drawing random cosmological populations of massive 
(AT > 5 x \0 u M Q h- 1 ) triaxial dark matter haloes at redshifts 
Z\ > 0.5 (cf. Appendix [C]). To investigate the impact of varia- 
tions of the cosmological parameters on the precise functional 
form of the correlation, we analysed cluster samples within both 
a WMAP1 and a WMAP7 cosmology. Furthermore, since the 
steepness of the inner core of dark matter density profiles is 
known to substantiall y change the s trong-lensing properties of 
galaxy clusters (see Bartelmann 2010), we additionally consid- 
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Fig. 7: Comparison between effective and median Einstein radii 
for a cosmological population of triaxial dark matter haloes with 
masses M > 5 x lO l4 M Q h~ l at redshifts z\ > 0.5, assuming 
a WMAP7 cosmology and adopting density profiles with inner 
slope a = 1.0. Each cross represents a pair of Einstein radii. 
The red dashed line indicates the bisector of the log (#E, m ed)- 
log (%eff) plane. The effective Einstein radii are systematically 
smaller, as is to be expected for ellipsoidal surface mass densi- 
ties. 



Fig. 8: Strong correlation between median Einstein radii %med 
and lensing cross sections cr 7 5 for a cosmological population 
of triaxial dark matter haloes (M>5xl0 14 M o r 1 ) at redshifts 
Z\ > 0.5, adopting a density profile with inner slope a - 1.0. The 
population was drawn assuming a WMAP7 cosmology. The red 
solid line indicates the best linear fit relation between log(cr 7 5 ) 
and log(#Emed)- The correlation breaks down for small Einstein 
radii (% med < 8"). 



ered the two extreme values of the inner slope, namely a = 1 .0 
and a = 1.5. 

Before we analyse the discussed correlation for the differ- 
ent settings, we begin with a short comparison of the alterna- 
tive definitions of Einstein radii. Figure [7 shows effective versus 
median Einstein radii for one specific population of dark mat- 
ter haloes, namely assuming a WMAP7 cosmology and adopt- 
ing triaxial density profiles with inner slope a = 1.0. We 
note that these plots qualitatively look the same for all other 
realisations (WMAP1 cosmology or steeper density profiles). 
Obviously, the agreement between the two definitions is mod- 
erate, and the absolute deviations increase with larger Einstein 
radii. Moreover, almost all points lie below the bisector, indicat- 
ing that effective Einstein radii are systematically smaller than 
median Einstein radii. However, this was to be expected, since 
both definitions only agree for axially symmetric lenses and start 
to differ (#E, e ff < #E,med) as soon as ellipticity is introduced. The 
deviation is especially pronounced when highly ellipsoidal sur- 
face mass density profiles are considered. 

In Fig. [8j we plot lensing cross sections cr 7 5 as a function 
of median Einstein radii #E,med for the same cluster sample as 



performed a least- squares fit for all data sets, assuming a linear 
relation of the form 



before. We can evidently reproduce the findings of Meneghetti 
et al.| ( |20lT] ). The correlation is especially tight for large Einstein 
radii, and thus also for large lensing cross sections. The scat- 
ter progressively increases as the size of the Einstein radii de- 
creases until the correlation finally breaks down for Einstein 
radii #E,med < 8". We recall that we considered elliptical sources 
with an equivalent radius of r s = 0.5" to compute lensing cross 
sections. In the regime where the correlation breaks down, the 
size of these sources is roughly of the same order as the size of 
the caustics, so that the efficiency to produce highly elongated 
images rapidly drops. Hence, as soon as the considered Einstein 
radii become too small, the correlation cannot be used anymore 
to infer lensing cross sections. Of course this breakdown also 
occurs in the case of effective Einstein radii. 

The correlation qualitatively looks the same for all other set- 
tings and using effective Einstein radii #E, e ff instead. Hence we 



log(cr 75 ) = alog(0 E ) + b . 



(29) 



Additionally, we computed Pearson's correlation coefficient r. 
The results of this procedure are summarized in Tab. [T] and sug- 
gest the following conclusions. First, we note that the coeffi- 
cients of the linear fits are not sensitive to the precise choice 
of the cosmological parameters. Furthermore, all Pearson cor- 
relation coefficients are close to unity, indicating that not only 
median, but also effective Einstein radii can be used to infer 
lensing cross sections in statistical studies. As may have been 
expected, the best-fittting linear relations between median and 
effective Einstein radii agree moderately. 

Next, steeper density profiles produce slightly larger Einstein 
radii, while the lensing cross sections increase substantially 
( |Oguri|[2004] ). This difference between steeper and flatter pro- 
files manifests itself in systematically larger intercepts b of the 
linear fits for steeper density profiles with inner slope a = 1.5. 
Furthermore, it explains why the Pearson correlation coeffi- 
cients indicate that the correlations are systematically stronger 
for a - 1.5. The distribution of Einstein radii is shifted towards 
higher values and, as already mentioned, the overall scatter of the 
correlation decreases as the Einstein radii increase. Additionally, 
fewer Einstein radii fall into the region where the correlation 
breaks down. Both effects can be verified in Fig. [9] which shows 
effective Einstein radii vs. lensing cross sections for the sam- 
ple assuming a WMAP7 cosmology and adopting an inner slope 
a - 1.5. Besides having different intercepts, we note that the 
slopes of the best linear fits are systematically flatter for a = 1.5. 
We conclude that steeper density profiles lead to flatter correla- 
tions between Einstein radii and lensing cross sections. We con- 
firmed this assumption with the following simple test. We gen- 
erated another random sample of triaxial dark matter haloes and 
placed all lenses at one single redshift, namely zi = 0.5. In that 
case, the geometrical configuration of the lens system is such 
that arcs occur farther away from the halo centres where the log- 
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Table 1 : Results of the least-squares fits for the relation between Einstein radii and lensing cross sections, assuming a correlation of 
the form log (cr 7 5) = a log (#e) + b. The last column shows the Pearson correlation coefficient r (see text for details). 



Cosmology 


Einstein radius 


a 


a 


b 


r 






1.0 


2.44 + 0.03 


-5.68 + 0.03 


0.96 


WMAP7 


#E,med 


1.5 


2.31+0.02 


-5.35 + 0.03 


0.98 




1.0 


2.40 + 0.04 


-5.54 + 0.04 


0.95 




#E,eff 


1.5 


2.29 + 0.02 


-5.26 + 0.03 


0.97 






1.0 


2.49 + 0.01 


-5.76 + 0.02 


0.96 


WMAP1 


#E,med 


1.5 


2.25 + 0.01 


-5.30 + 0.01 


0.98 




1.0 


2.46 + 0.02 


-5.63 + 0.02 


0.95 




#E,eff 


1.5 


2.26 + 0.01 


-5.24 + 0.01 


0.98 
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Fig. 9: Strong correlation between effective Einstein radii #E, e ff 
and lensing cross sections cr 15 for a cosmological population 
of triaxial dark matter haloes (m>5x 1O 14 M /z _1 ) at redshifts 
Z\ > 0.5, adopting a density profile with inner slope a = 1.5. The 
population was drawn assuming a WMAP7 cosmology. The red 
solid line indicates the best linear fit relation between log(<T7.5) 
and logoff). 



arithmic slope of the density profile is steeper. Consequently, we 
found a substantially flatter correlation for that specific sample. 



numerical clusters 




#E,med [arcsec] 

Fig. 10: Best linear fits of the correlation between median 
Einstein radii #E,med and lensing cross sections <x 7 .5 for numer- 
ical and analytic clusters. The blue solid line indicates the best 
linear fit for numerical clusters (cf. Eq. ([28])). The red dashed and 
the green dotted line indicate the best linear fits for cosmological 
populations of triaxial dark matter haloes (cf . Table [I]) with in- 
ner profile slopes a = 1.0 and a - 1.5, respectively, assuming a 
WMAP1 cosmology. The lensing cross sections of triaxial dark 
matter haloes are systematically larger in the range #E,med > 8". 



4.4. Comparison between numerical and analytical clusters 



We n ow compare our results to the findings of Meneg hetti et al. 
( 2011) in detail. As indicated by Fig.fTOl the lensing cross sec- 
tions of triaxial haloes as a function of median Einstein radii 
are systematically larger than those of the numerical clusters 
throughout the entire range where the correlation between both 
quantities holds (#E, m ed ^ 8"). This discrepancy even increases 
as the gravitational lenses become stronger. Since there is no ob- 
vious explanation for the above deviations, we briefly discuss 
the main differences between the two cluster samples. First, the 
MareNostrum simulation contains baryons, which are known to 
enhance the lensing cross sections of galaxy clusters ([Puchwein 
eTaL][2005| |Wambsganss"eraLl[2008l |Rozo et al.||2008[ |Mead 
et al. 2010| . Second, we describe dark matter haloes as ellip- 
soids with constant axis ratios, although numerical simulations 
indicat e that these ratios slightly decrease towards the halo cen- 
tres ( [Jing & Suto 2002). Hence we slightly underestimate the 
ellipticity in the inner part of the density profile. Third, numer- 
ically simulated clusters naturally contain substructure and in- 
corporate the effect of cluster mergers, whereas we only consid- 



ered isolated dark matter haloes. All these previous aspects, i.e. 
asymmetry, substructure, and cluster mergers, are known to sub- 
stantially increase the strong-lensing efficiency of galaxy clus- 
ters flBartelmann et al.|1995| |Torri et al.|2004| |Fedeli et al.|2 006 ; 



Meneghetti et al. 2007). 



Given these differences, the lensing cross sections of numer- 
ical clusters should be significantly larger. We speculate that the 
inclusion of substructure and cluster mergers is the main differ- 
ence between the cluster samples. If this is indeed the case, Fig. 



10 suggests that these cluster properties become increasingly im- 



portant for stronger gravitational lenses. Substructure leads to 
irregular critical curves with some extreme outliers. These, how- 
ever, should be compensated for by the computation of median 
Einstein radii. Hence, we expect that substructure, on average, 
increases the median Einstein radii only slightly, whereas lens- 
ing cross sections are likely significantly enhanced. Evidently, 
the comparison of the correlations in Fig. [10| indicates exactly 
the opposite trend, implying that substructure alone cannot ex- 
plain the deviations. On the other hand, the evolution of lensing 
cross sections as a function of median Einstei n rad ii during clus- 
ter mergers is far from obvious. In Section 5.3 we show that 
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these eve nts might indeed explain the discrepancy between the 
results of Meneghetti et al. (201 1) and our findings. 



5. The evolution of strong-lensing properties during 
cluster mergers 

5.1. Evolution of tangential critical curves 

We now discuss the evolution of a tangential critical curve dur- 
ing a cluster merger by means of the following simplified toy 
model (cf. |Torri et al.|2 004 ). We considered two massive triaxial 
haloes merging at the fixed redshift z\ = 0.5. The source plane 
was placed at z s = 2.0. The main halo (M main = 1 x 10 15 M Q h~ 1 } 
was assumed to rest at the coordinate origin, while the second 
halo (M sub = 2.5 x 10 14 M Q h~ l ) was placed at an initial distance 

d ~ 1 Mpc h~ l from the main halo. Since the precise prop- 
erties (concentration, axis ratios, orientation) of both haloes do 
not matter for the following discussion, we simply state that their 
shape and orientation were chosen such that both lensing po- 
tentials are considerably elliptical. We simulated the merger by 
successively reducing the separation between both haloes in a 
series of discrete steps until their centres finally overlapped. At 
each step, we computed the resulting tangential critical curves. 



Figure [TT] illustrates the evolution of the tangential critical 
curves for a selection of intermediate steps of the simulated 
merger. As can be seen, we oriented the main halo such that 
the semi-major axis of its elliptical surface mass density pro- 
file was aligned with the direction of motion, whereas the ellip- 
tical projected mass profile of the second halo was rotated by 
90 degrees. In the beginning (top left panel), the separation be- 
tween the haloes is large so that both have their own isolated 
and almost unperturbed tangential critical curves. As soon as the 
haloes approach, the shear in the region between them grows 
and stretches both tangential critical curves along the direction 
of motion until they finally merge to build one large, highly elon- 
gated structure (d^0.32 Mpc /T 1 ). Given the tight correlation 
between Einstein radii and lensing cross sections, we expect that 
these configurations are particularly efficient in producing thin 
arcs (see next section). As the smaller halo continues to approach 
the main halo, the highly elongated tangential critical curve typi- 
cally starts to shrink along the direction of motion. Finally, how- 
ever, the two mass profiles start to overlap significantly so that 
the convergence of the inner region grows substantially. As a 
consequence, the tangential critical curve again starts to expand 
almost isotropically. 

This interpretation of the evolution of critical curves during 
cluster mergers was first given by |Torri et aT] ( |2004| ), and we 
refer the reader to their work for more details. It goes without 
saying that there are more sophisticated studies of cluster merg- 
ers that take virialization and other important physical processes 
properly into account. However, these models clearly go beyond 
the scope of our fast, semi-analytic approach to investigate the 
statistical strong-lensing properties of large cosmological popu- 
lations of triaxial dark matter haloes. 



5.2. Evolution of Einstein radii and lensing cross sections 



clusters, but also for the overall statistical strong-lensing effi- 
ciency of the cosmic cluster population. Motivated by these re- 
sults, we performed a series of tests using the above toy model 
to study (1) whether we are able to reproduce the increase of the 
strong-lensing efficiency during mergers by modelling clusters 
with triaxial density profiles, (2) which parameters dominate the 
strength of this effect, and (3) whether the correlation between 
Einstein radii and lensing cross sections still holds during these 
events. 

Our most important results can be summarized as fol- 
lows. Naturally, a notable enhancement of the strong-lensing 
efficiency can only be observed during major mergers (i.e. 
^sub/^main » 0.05). If the mass ratio M su b/M ma i n is too low, 
the elongation of the tangential critical curve enclosing the main 
halo is negligible. Conversely, the effect is particularly strong for 
mergers of clusters with comparable mass. Furthermore, the rel- 
ative orientation of the ellipsoids plays a dominant role. Figure 
12] sketches four particular relative orientations of merging tri- 



Torri et aL] ([2004) numerically simulated cluster mergers and 



found that these events are capable of boosting the lensing cross 
sections for giant arcs by about one order of magnitude. In a sub- 
sequent work, Fede li et aL] ([2006 ) showed that these events are 
not only important for the lensing cross sections of individual 



axial dark matter haloes. The blue ellipsoid in the middle rep- 
resents the projected mass density of the main halo, while the 
red ellipsoids display the surface mass densities of the infalling 
subhaloes. As discussed above, the enhancement of the strong- 
lensing efficiency can be observed shortly before the two ini- 
tially separated tangential critical curves merge to form one 
large, extremely elongated structure. This effect is especially 
pronounced when the direction of motion and the semi-major 
axis of the main halo's surface mass density profile are aligned 
(mergers (B) and (D) in Figure [12]). The most extended tangen- 
tial critical curves form in situation (B), where the semi-major 
axes of both projected mass profiles are aligned and addition- 
ally point along the direction of motion. Then, the elongation 
of the two separated tangential critical curves sets in relatively 
early and the diameter of the merged critical curve is particu- 
larly large. We emphasize that we can reproduce an increase of 
the lensing cross section by factors of ~ 2 for typical mass ratios 
^sub /^main ~ 0.25 if, and only if, this distinguished merger con- 
figuration is chosen and the surface mass density profiles of both 
haloes are appreciably elliptical (see below). In all other cases, 
the enhancement caused by mergers is notably smaller. Clearly, 
we cannot reproduce the order-of-magnitude increase found by 
|Torri et aL] (|2004). Particularly in the least favourable configura- 
tion (C) there is almost no enhancement at all. 

To quantify these general remarks, we exemplarily consider 
two specific cluster mergers in detail. In both cases, we simu- 
lated triaxial haloes of masses M ma i n = 1 x 10 15 M Q h~ l and 
M sub = 2.5 x 10 14 M Q h~ l , with axis ratios a/c = 0.4 and 
b/c = 0.6, mean concentration c Q « 1.5 (cf. Eq. fTO])), inner 
profile slope a = 1.0 and projected along the / -direction. 

We first consider a merger choosing orientation (D). Figure 
[T3l shows the evolution of different definitions of Einstein radii 
and the lensing cross section as a function of the distance be- 
tween the haloes' centres. In both diagrams, the dashed vertical 
line indicates the moment when the initially separated tangential 
critical curves merge. The blue solid and the red dashed lines in 
the left panel indicate the effective and median Einstein radii, re- 
spectively, computed for the tangential critical curve enclosing 
the main halo only. These definitions take into account that both 
haloes initially have their own independent critical structures and 
are thus distinguishable systems from the lensing point of view. 
In contrast, the green dotted line indicates the median Einstein 
radius of all tangential critical points relative to the centre of 
the main halo. Here, "all" means that for the larger distances, 
where both clusters have their own critical curves, also the criti- 
cal points enclosing the merging subhalo are taken into account 
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Fig. 1 1 : Evolution of the tangential critical curves during a merger of two massive 
(M main = 1 x 10 15 M e h~\ M sub = 2.5 x 10 14 M Q h~ l ) triaxial dark matter haloes. Starting with an initial separation of 

d = 0.57 Mpc h~ l in the top left panel, the distance between the two haloes is successively reduced until their centres 
finally overlap in the bottom right panel. 



(A) 



(D) 




(B) 



(C) 



Fig. 12: Schematic sketch of four particular orientations of merg- 
ing triaxial dark matter haloes. The blue ellipsoid in the middle 
represents the projected mass density profile of the main halo. 
The red ellipsoids visualize the surface mass density profiles of 
the inf ailing subhaloes, which have different relative orientations 
with respect to the main halo. The black arrows indicate the di- 
rections of motion of the inf ailing structures. 



for computing the median Einstein radius. This definition closely 
resembles that used by Meneghett i et al. ( |2011| , who measured 
the distances of all tangential critical points with respect to the 
highest mass peak. We note that the effective Einstein radius 
seems to capture the isotropic expansion of the tangential criti- 
cal curve as the density profiles of the two haloes start to overlap 
significantly and the convergence of the central regions grows. 
However, it does not indicate a substantial boost of the strong- 
lensing efficiency when the critical curves merge. This effect is 
better captured by median Einstein radii, which on the other 
hand do not reflect the isotropic growth of the tangential criti- 



cal curve for short distances. The right panel of Fig. 13] shows 
the evolution of the total lensing cross section cr 7 5 of the two 
clusters. Obviously, the strong-lensing efficiency grows substan- 
tially shortly before the critical curves merge. For that specific 
orientation, we observe an increase of ~ 170% in comparison to 
the initial value. After that peak, the lensing cross section passes 
through a minimum - reflecting the fact that the elongation of the 
tangential critical curve decreases and the corresponding caustic 
becomes less cuspy - before it finally increases by ~ 200% as a 
consequence of the overlapping mass profiles. 

For comparison, Fig. [14] shows the same diagrams for orien- 
tation (A), which corresponds to a direction of motion perpen- 
dicular to the major axis of the main halo's surface mass density 
profile. Evidently, the evolution of the strong-lensing properties 
differs substantially from the previous merger. The elongation of 
the critical curve surrounding the central halo sets in later and 
the resulting Einstein radii are smaller. Moreover, the separa- 
tion between the two clusters needs to be smaller for the critical 
curves to merge. The first peak of the lensing cross section is 
significantly smaller and coincides with the increase caused by 
the growing mass concentration in the central regions. 



5.3. Correlation between Einstein radii and lensing cross 
sections 

Our main motivation for studying the evolution of Einstein radii 
and lensing cross sections is to investigate whether or not the 
correlation between the two quantities still holds during merg- 
ers. If the relation is still valid, we can later simply compute all 
Einstein radii (#e > 8") on the sky and infer their lensing cross 
sections by means of the correlation instead of actually com- 
puting them with one of the discussed methods (cf. Section [5]). 
This approach would save much computing time. This final ap- 
plication may also justify why we plotted different definitions of 
Einstein radii in the left panels of Fig. 13 and 14 since the first 
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Fig. 13: Figure showing the evolution of different definitions of Einstein radii (see text) and the lensing cross section during a 
simulated cluster merger, assuming configuration (B) in Figure 12 All quantities are plotted as a function of the distance between 
both halo centres. The dashed vertical line indicates the moment when the tangential critical curves merge. 
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Fig. 14: Evolution of different definitions of Einstein radii (see text) and the lensing cross section during a simulated cluster merger, 
assuming configuration (A) of Figure 12 All quantities are plotted as a function of the distance between the two halo centres. The 
dashed vertical line indicates the moment when the tangential critical curves merge. 



two of them only consider the tangential critical curve enclos- 
ing the main halo and are hence clearly more appropriate for our 
purpose. 

The following analysis focuses only on the regions where 
both critical curves have already merged (i.e. left of the dashed 
vertical lines in Figs. [13] and [14]). First, we observe that none of 
the definitions of Einstein radii is able to preserve the tight corre- 
lation with lensing cross sections. In any case, there are regions 
where the Einstein radii decrease while the lensing cross sections 
increase, and vice versa. The correlation seems to break down 
completely if median Einstein radii are used. In comparison, ef- 
fective Einstein radii trace the evolution of lensing cross sections 
better. However, by our definition of measuring only the area that 
is enclosed by the tangential critical curve surrounding the main 
halo, they do not capture the enhancement of the strong-lensing 
efficiency shortly before the tangential critical curves merge. We 
observe in passing that the definition of effective Einstein radii is 
more natural for cluster mergers, because it does not require an 
arbitrarily chosen fiducial point. 



In Sect. |5.1[ we assumed a constant lens redshift of z\ = 0.5 
for our simplified models of cluster mergers. For reference, we 
therefore again generated a sample of 300 random triaxial haloes 
(mass range 5xl0 14 -2xl0 15 M Q h~ l , inner profile slope a = 1.0) 
placed at zi = 0.5 and computed the best linear fits for the cor- 
relation between the two different definitions of Einstein radii 
and lensing cross sections. Figure [T5| plots these fits, where the 
left and the right diagrams correspond to median and effective 
Einstein radii, respectively. In addition, we simulated all four 
merger configurations of Fig. 12 using the same dark matter 
haloes as in the previous section and computed their Einstein 
radii and lensing cross sections at each discrete simulation step. 
The blue crosses in the left and right diagrams indicate the re- 
sults of all those discrete steps that correspond to distances with 
merged tangential critical curves. Evidently, for median Einstein 
radii, the pairs (0E, m ed> o~i.s) systematically lie below the corre- 
lation found for isolated clusters. This can easily be understood 
from Figs. [13] and [14] Starting at d = 0.0, the median Einstein 
radii increase as a function of distance, whereas the lensing cross 
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sections decrease, so there is indeed an anti-correlation. This ob- 
servation provides a possible explanation for the discrepancy be- 
tween the linear relation found by Meneghett i et al.| ([201 1 ) and 
our results (cf. Section |4~4| ). We expect that a certain fraction of 
Meneghetti and coworkers' numerically simulated clusters are 
actively merging so that their contribution flattens the results of 
the least- squares fit. If that line of reasoning is correct, we expect 
that the correlations found for numerical clusters predict smaller 
lensing cross sections as a function of Einstein radii. It would be 
an interesting task to verify this assumption by studying selected 
cluster samples in numerical simulations. The right-hand side di- 
agram of Fig. [15] shows that the pairs of effective Einstein radii 
and lensing cross sections during cluster mergers agree signifi- 
cantly better with the corresponding linear fit found for isolated 
dark matter haloes. Hence we conclude that effective Einstein 
radii capture the evolution of strong-lensing properties during 
cluster mergers better and the correlation can still be used to ap- 
proximate lensing cross sections in statistical lensing studies. 



6. Impact of cluster mergers on the statistics of the 
largest Einstein radii 

6.1. A semi-analytic approach to sample large cosmological 
halo populations that incorporates cluster mergers 

While large cosmological Af-body simulations are probably the 
most realistic framework for studying the impact of mergers on 
the statistical lensing properties of selected cluster samples, they 
are computationally costly. As an alternative, we now introduce 
a fast, semi-analytic method that allows us to study the statisti- 
cal impact of mergers within a fraction of the computing time 
required for Af-body simulations. We follow this approach be- 
cause we aim to be able to repeat our calculations for large 
cosmic volumes, varying a multitude of boundary conditions 
(i.e. different properties of the triaxial density profile, cosmo- 
logical parameters, mass functions, etc.), within a reasonably 
short time. This ability is important for our follow-up work, 
where we plan to employ methods of extreme value statistics 
to investigate the probability distribution of th e largest Einstein 
radius expected within a ACDM cosmology ( [Waizmann et al. 
|2012| ). We require our new approach to be fully self-consistent 
within the framework of the extended Press-Schechter theory 
( Lacey & Cole|[T993] ), which means that in particular the pre- 
dicted merger rates and t he Press-Schechter mass function (cf. 
Press & Schechter 1974| ) must be reproduced accurately. This 



way, our results should be comparable with those derived from 
numerical simulations and real observations. 

Let us assume that we intend to analyse the strong-lensing 
properties of a representative all- sky realization of dark mat- 
ter haloes in the mass and redshift range [M m i n ,M max ] and 
[Zmm, ZmaxL respectively. To do this, we sequentially perform the 
following steps: 

1 . Initial halo population: We use a common Monte Carlo ap- 
proach (cf . |Oguri & Blandford| ( |2009| ) and Appendix [C) to 
populate the considered cosmic volume with an initial sam- 
ple of dark matter haloes at the present time (z = 0). We sub- 
divide the mass range into logarithmically equidistant bins 
of typical size A(log(M)) = 0.02 and calculate the mean 
expected number N of haloes in each bin using the Press- 
Schechter mass function dn(M, z = 0)/dM, 



Here, V(z m in,£max) denotes the comoving volume of the 
spherical shell between the minimum and the maximum 
redshift (with respect to an observer at the coordinate ori- 
gin). We then generate a random integer number N from 
the Poisson distribution with mean N and sample N haloes 
with logarithmically uniformly distributed masses in the cor- 
responding bin. Since strongly lensing clusters are rare ob- 
jects that are generally separated by large distances, we can 
neglect any kind of large scale structure correlation and sim- 
ply uniformly distribute these N sampled haloes in the con- 
sidered cosmic volume. To this end, we assign each halo a 
random position on the sky, generate a uniformly distributed 
random number x in the range [0, 1) and then compute the 
comoving distance r with respect to the observer, 



[*r 3 (z max ) + (1 



x)r 3 (Zn 



.)] 



1/3 



(31) 



where r(z m m) and r(z max ) denote the comoving distance to 
the minimum and maximum redshift, respectively. This for- 
mula takes the varying volume (as a function of radius) of 
spherical shells into account. Given the comoving radius r, 
we can easily infer each halo's redshift z bs with respect to 
the observer. 

Reverse time evolution: We evolve the initial halo popula- 
tion backwards in time. More precisely, we adop t the effi- 
cient Monte Carlo approach proposed by Zhang et al. ( 2008 



see 'method B') to simulate a representative merger tree of 
each individual halo up to its previously de termined observa- 
tion redshift z bs- Since Zhang et al. (2008) provide a nicely 
written step-by-step description of 'method B', we refer the 
reader to their work for details. Generally, merger tree algo- 
rithms generate representative formation histories by evolv- 
ing the initial halo mass backwards in time taking discrete 
time steps (typical step size Az ~ 0.02). At the first time step, 
the halo is split into sm aller progenitors by means of the con- 
ditional mass function ( [Lacey & C ole 1993). At the next time 
step, these new progenitors are again split into yet smaller 
progenitors using the same recipe, and so forth. Applying 
that scheme to all arising progenitors, this finally yields the 
full, tree-like formation history (i.e. the merger tree) of the 
initial halo. Given our predefined mass range, however, we 
discard all progenitors with masses below M^. We note that 
we verified that our implementation accurately reproduces 
the Press-Schechter mass function as well as the theoreti- 
cally predicted merger rates at any look-back time. 
3. Kinematics of merger trees: To follow the trajectories of 
all arising progenitors in the considered cosmic volume, we 
additionally need to describe the kinematics of the merger 
trees. For that purpose, we adopt the following simplistic ap- 
proach, which was first introduced by Fedeli & Bartelmann 
( 2007| ). Each time a merger between two haloes of masses 
Mi and M 2 occurs, we estimate its duration using the dy- 
namical time scale 



yn 



/ Ovir,l + 7Vir,2) 

' G (M\ + Mi) 



(32) 



N = 



dn(M,z = 0) 
dM 



AM X Vfonin, Zmax) • 



(30) 



where r v i r? i and r v i r? 2 denote the virial radii of both haloes. 
These dynamical time scales are typically of the order of 
a several hundred Myr, which agrees well with merger du- 
rations measured in numerical simulations (see |Torri et aL] 
2004 for instance). Assuming a uniform linear motion, we 
compute the relative velocity of the two haloes, 

^vir,l + r \ir,2 

Vrel = , (33) 

* dyn 
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Fig. 15: Comparison of the correlation between Einstein radii and lensing cross sections for a sample of isolated triaxial dark matter 
haloes to the same correlation for a sample of merging triaxial dark matter haloes. The red solid line in the left (right) diagram 
indicates the best linear fit relation between median (effective) Einstein radii and lensing cross sections found for a sample of 
triaxial dark matter haloes at redshift zi = 0.5. The blue crosses represent selected pairs of Einstein radii and lensing cross sections 
during cluster mergers (see text for details). 



and finally sample a random direction of motion. We apply 
this scheme to all mergers. This way, we can easily com- 
pute the spatial positions of all progenitors of the considered 
merger tree at the observation redshift z bs- We emphasize 
that our model of cluster mergers is clearly simplistic and 
neglects certain physical details. For instance, we implic- 
itly assumed central cluster collisions although the centres of 
merging clusters can t ypically be sep arated by up to several 
hundred kpc (see e.g. [Sarazin 2002| for a semi-analytic ap- 
proximation). Finite impact parameters may certainly influ- 
ence the detailed evolution of the strong-lensing properties 
during individual cluster mergers, and it would be interesting 
to study refinements of our simplified model in future works. 
However, we do not expect that these simplifications have 
a significant impact on the cumulative lensing efficiency of 
cosmological samples of merging clusters. 

Summarizing the above steps, we first populate the considered 
cosmic volume with an initial halo sample at redshift z = 0, 
evolve that sample backwards in time using merger trees and a 
simplistic model to describe the kinematics, and finally project 
the resulting halo configurations onto the observer's past null 
cone. This procedure is illustrated in Fig. [T6] 

We conclude this section with a final remark on the spe- 
cific choice of the mass function. Our algorithms are based on 
the Press-Schechter mass function ( [Press & "S chechter| |1974| ) 
throughout, although improved variants of the mass function 



with better accuracy have been proposed (see Jenkins et al. 
ISheth & Tormenl|j502l |Warren et al.1[2006l [Tinker et al 



2001 



2008 



for instance). This choice is motivated by the fact that the Press- 
Schechter mass function is based on the theory of spherical col- 
lapse and only within this framework exact analytic expressions 
for the conditional mass function - which is the fundamental 
quantity for the computation of merger trees - can be derived. 
In contrast, most improved variants are empirical fits to numer- 
ical simulations and hence there is no theoretical framework to 
derive the corresponding expressions for the conditional mass 
function. Since we aim to compare the statistics of Einstein radii 
of samples of isolated haloes to those of halo samples that in- 



corporate cluster mergers in a fully self-consistent way, we are 
restricted to the mass functions based on the spherical (or ellip- 
soidal) collapse. The moderate deviations of the Press-Schechter 
mass function are negligible for our principal conclusions in this 
work. In contrast, if we were to perform a detailed comparison 
between our theoretical predictions and observations, the pre- 
cise choice of the mass function would play an important role. 
This is mainly because the strongest gravitational lenses stem 
from the exponentially suppressed high-mass tail of the mass 
function, which is particularly sensitive to the exact calibration. 
Because the Press-Schechter mass function is known to under- 
estimate the abundance of haloes in the high-mass tail, our ap- 
proach also slightly underestimates the number of extreme lens- 
ing events. We refer the interested reader to our follow-up work 
( Waizmann et al.|2012| ), where we provide a detailed discussion 
of the impact of different mass functions on the statistics of the 
single largest Einstein radius. 



6.2. Computation of Einstein radii 

The previous step populates the observer's past null cone with 
(merging) dark matter haloes. We now describe an efficient 
method for computing the distribution of Einstein radii of that 
halo population. 

First, we sort the masses in descending order and begin with 
the computation of the Einstein radius of the most massive halo 
Mi . To this end, we scan the region around that halo for neigh- 
bouring subhaloes M t whose distance d is shorter than the sum 
of both virial radii, d < (r v i r j + r v i r?i ). If any surrounding sub- 
haloes satisfy this condition, they are taken into account for the 
computation of deflection angles. All other haloes that are far- 
ther away do not influence the strong-lensing properties of the 
considered massive halo M\ and can safely be ignored. Next, 
we assign a random triaxial density profile and a random ori- 
entation with respect to the observer to each relevant halo. This 
enables us to determine the tangential critical curve that encloses 
the massive halo M\. To do this, we place halo M\ at the coor- 
dinate centre of a Cartesian grid whose side length is chosen to 
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Fig. 16: Sketch illustrating our new algorithm for generating re- 
alistic cosmological halo populations that incorporates clusters 
mergers. An initial halo population at redshift z - is evolved 
backwards in time using extended Press-Schechter merger trees 
and is finally projected onto the observer's past null cone. 



be sufficiently long to contain all relevant subhaloes. Instead of 
computing the deflection angle map of the complete field and 
determining the tangential critical curve afterwards, we simply 
need to detect the first tangential critical point left of the coor- 
dinate origin and employ a standard friend-of-friend algorithm 
to identify the entire tangential critical curve. This way, we only 
have to compute deflection angles in small stripes surrounding 
the critical curve, which renders the algorithm very fast. If we 
aimed to compute the lensing cross section of the considered 
configuration instead, we would not be able to use this efficient 
method but instead would have to calculate the complete de- 
flection angle map. This important difference illustrates why the 
determination of Einstein radii is computationally far less ex- 
pensive and why it is of great advantage to infer lensing cross 
sections from Einstein radii using the discussed correlation be- 
tween both quantities. Finally, we detect all haloes which are en- 
closed by the computed tangential critical curve and mark them 
as haloes which do not have an independent Einstein radius. We 
repeat the above procedure for all remaining, unmarked haloes 
of the sorted mass list so that we end up with a complete cata- 
logue of Einstein radii. 

6.3. Comparison of the statistics of Einstein radii with and 
without duster mergers and inferred optical depths 

We can now quantify the statistical impact of cluster mergers 
on the strongest gravitational lenses (#E, e ff ^ 10") in our uni- 
verse. To this end, we generated a first sample of single dark 
matter haloes using the conventional Monte Carlo approach (see 
Appendix [C]) and a second sample using our newly dev eloped 
technique to incorporate cluster mergers (see Section 6.1 ). Since 
the largest observed Einstein radii, which may challenge the 
standard cosmo logical model, were all found at redshifts wel l 
below unity (cf. [Broadhurst & Barkana|2008[|Zitrin et al.|20lT 



we focused our analysis on the redshift range z e [0, 1]. In an- 
ticipation of the following results, we note that the masses of all 



single haloes with effective Einstein radii above 10" are larger 
than 10 14 M@/h. Recalling that the mass ratio M su b/M ma i n of 
two merging clusters at least needs to exceed 5% to notably 
boost the strong-lensing efficiency of the main halo, we could 
therefore safely ignore all haloes with masses below M m \ n = 
5 x 10 12 MQ/h. In particular, we pruned all branches (i.e. pro- 
genitors) of the merger trees whose masses dropped below that 
threshold. In this section, we modelled all haloes by means of 
triaxial density profiles with an inner slope a = 1, adopted a 
WMAP7 cosmology and assumed that the source plane is placed 
at redshift z s = 2.0. 

Figure 17 and Tab. [2] clearly demonstrate that cluster merg- 
ers do have a significant impact on the distribution of the largest 
Einstein radii. While we find 4132 Einstein radii larger than 
10" in the single halo run, that number increases by 36% up to 
5622 haloes if mergers are taken into account. The effect is even 
more significant if we only consider the number of systems with 
Einstein radii above 20", in which case we find an increase of 
74%. Another interesting aspect to be analysed is the frequency 
of cluster mergers. To do this, we classified observed systems as 
being actively merging if more than one halo is enclosed by the 
tangential critical curve. Clearly, this definition does not cover 
those cases where merging clusters are already close to each 
other, but still have their own, highly elongated critical curves. 
Nevertheless, our wantedchoice should be sufficiently accurate 
for a qualitative estimate. We find that 35% of the systems with 
Einstein radii larger than 10" are actively merging. This num- 
ber increases to 55% if we only consider those clusters with 
Einstein radii larger than 20". Given these results and the above 
number counts of Einstein radii, we can conclude that (1) clus- 
ter mergers are an important mechanism to increase the statisti- 
cal lensing efficiency of cosmological halo populations and that 
(2) these events become increasingly dominant for particularly 
strong gravitational lenses. 

Additional evidence in favour of these conclusions is pro- 
vided by considering the very largest lenses in detail. In the sim- 
ulation including mergers, eight out of the ten largest Einstein 
radii stem from actively merging systems. Moreover, the largest 
Einstein radius in the single halo run has the size #E, e ff = 38.5", 
while we find a notably larger maximum of #E, e ff = 50.8" in 
the run with mergers. Needless to say, these maxima are subject 
to statistical fluctuations. Since the considered halo populations 
were randomly drawn, it may well be that for other realizations 
the single halo maximum is larger than the one found with merg- 
ing haloes. However, we can safely conclude that cluster mergers 
need to be taken into account in semi-analytic, statistical stud- 
ies that aim to challenge the standard cosmological model on 
the basis of the largest observed Einstein radii. The left panel of 
Fig. 18 shows the critical curve of the largest observed Einstein 
radius, which belongs to an extraordinary system at redshift 
z = 0.43. Starting with an initial mass of M = 1.8 x 10 15 M Q /h 
at redshift z = 0, the halo was split into 20 progenitors with 
masses above the minimum mass threshold M m i n before it fi- 
nally reached the observer's past null cone. Six of these 20 pro- 
genitors are enclosed in the critical curve shown. The most mas- 
sive halo in the centre has a mass of M = 1.1 x 10 15 Mq/Ji. 
The masses of the five surrounding subhaloes lie in the range 
6.7 x 10 12 - 1.5 x 10 14 M o //l To further demonstrate the im- 
portance of mergers, we removed these subhaloes and computed 
the Einstein radius of the main halo alone, finding that it drops 
by ~ 20" to now only #E,eff = 30.9". Moreover, we can use 
that system to verify our choice of the minimum mass thresh- 
old M min = 5.0 x 10 12 M©//z. The smallest enclosed subhalo 
has a mass of M = 6.7 x 10 12 M©//i, which is close to M m i n . If 
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we ignore that halo for the computation of the tangential criti- 
cal curve, the Einstein radius drops only slightly by ~ 1 .4% to 
#E,eff = 50. 1", confirming that our mass cut-off only leads to 
negligible errors. There are mainly two reasons for introducing 
such a mass cut-off in general: (1) we can discard a consider- 
able amount of small-mass branches of the merger trees and (2) 
we can limit the number of (irrelevant) subhaloes that need to 
be taken into account for the computation of the critical curves. 
Both aspects substantially improve the performance of our algo- 
rithm. 

The framework of hierarchical structure formation suggests 
that relatively young clusters at intermediate and high redshifts 
(z > 0.5) are dynamically more active than older systems at low 
redshifts. Therefore, it is to be expected that mergers mainly in- 
crease the number of strong lenses in the upper redshift range of 
our simulation volume. Although Fig. 19 confirms this expecta- 
tion, the net shift of the overall population of strong-lensing clus- 
ters to higher redshifts due to mergers is only moderate. Instead, 
the redshift distribution of the strongest lenses is predominantly 
determined by the chosen lensing geometry with sources at red- 
shift Zs = 2.0. 

Finally, we used the correlation between Einstein radii and 
lensing cross sections to infer the impact of cluster mergers 
on the statistics of giant gravitational arcs from the distribu- 
tions of Einstein radii just computed. Since we concentrated on 
the strongest gravitational lenses only, we are well inside the 
regime where the correlation between both quantities is partic- 
ularly tight and hence our following estimate should be reason- 
ably accurate. First, we randomly picked ~ 400 haloes of our 
catalogue with Einstein radii above 10" and computed their lens- 
ing cross sections cr 7 5 . We then used these results to calibrate the 
best linear fit relation for the correlation: 

log(<r 7 . 5 ) = (2.12 ± 0.06) logoff) - (5.24 ± 0.07) . (34) 

The cumulative lensing efficiency of cluster samples is usually 
characterized by their optical depth r for giant gravitational arcs, 
which is given by the sum of the individual lensing cross sections 
cr 15 i divided by the size of the entire source sphere, 



T = 



(4nD 2 s y l 



(35) 



where D s denotes the angular diameter distance to the source 
plane. Assuming that the individual lensing cross sections do 
not overlap in the source plane, the optical depth corresponds to 
the probability that an arbitrarily placed source with the speci- 
fied characteristics produces a gravitational arc with a length-to- 
width ratio higher than 7.5. We find that cluster mergers increase 
the optical depth of all haloes with Einstein radii above 10" by 
roughly 45%. Furthermore, as was to be expected, the impact on 
those systems with Einstein radii above 20" is even stronger: we 
find that mergers statistically increase the number of giant arcs 
produced by these particularly strong lenses by roughly 85%. 

7. Conclusions 

Our main results can be briefly summarized as follows. As a first 
preparatory step, we carried out a detailed comparison between 
three alternative methods to compute lensing cross sections. We 
sh owed that the fast, s emi-analytic method originally proposed 
by |Fedeli et aL| p006| ) systematically deviates from ray-tracing 
algorithms due to a differing definition of the length-to-width 
ratio of lensed images. However, this discrepancy can easily be 



Table 2: Comparison of the strong-lensing statistics of two repre- 
sentative samples of dark matter haloes (M > 5x 10 12 Mq/Ii, z e 
[0, 1]) including and excluding cluster mergers, adopting a 
WMAP7 cosmology and a constant source redshift z s = 2.0. We 
show the maximum effective Einstein radius max(#E,eff ) includ- 
ing the observation redshift z bs of the corresponding lens sys- 
tem, the number of Einstein radii above a certain threshold X, 
N (#E, e ff ^ X), and the optical depth (for gravitational arcs, see 
Eq. ( [35] )) of all haloes with Einstein radii above the threshold X, 
T(0E,eff > X). 





single haloes 


merging haloes 


max(% eff ) 


38.5" (z obs = 0.21) 


50.8" (z obs = 0.43) 


Af(%eff> 10") 


4132 


5622 (+36%) 


W(%efF>20") 


133 


231 (+74%) 


T(%eff > 10") 


2.9 x 10" 7 


4.2 x 10" 7 (+45%) 


T(% eff > 20") 


3.3 x 10" 8 


6.1 x 10" 8 (+85%) 



with mergers (N=5622) 
active mergers (N=1968) i i 
without mergers (N=4132) Y////A 




0.4 0.6 
redshift z 

Fig. 19: (Normalized) redshift distribution of the strong gravi- 
tational lenses with effective Einstein radii larger than 10" in- 
cluding (blue solid histogram) and excluding (black patterned 
histogram) cluster mergers. In addition, the red stacks indicate 
the fraction of actively merging systems of all strong gravita- 
tional lenses represented by the blue solid histogram, i.e. includ- 
ing cluster mergers. Systems are classified as actively merging 
if more than one halo is enclosed by the tangential critical curve 
(see text). 



resolved by a trivial correction, after which the semi-analytic 
method yields equally reliable results as the well-probed ray- 
tracing codes. Hence, our findings confirm that the semi-analytic 
method represents an extremely valuable tool for computation- 
ally demanding problems. Moreover, our analysis clearly indi- 
cates that the precise definition of the length-to-width ratio of 
lensed images has a significant impact on the final lensing cross 
section. Thus, special care must be taken to adopt consistent 
definitions, in particular when comparing theoretical predictions 
with real observations in the context of the arc statistics problem. 

Secondly, we investigated the correlation between Einstein 
radii and lensing cros s sections that was found recently by 
Meneghetti et al. ( 201 1) for a selected sample of numerical clus- 
ters. Considering single haloes only, we showed that we can re- 
produce the tight correlation, using both median and effective 
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Fig. 17: Comparison of the distributions of the largest Einstein radii including (blue patterned histogram) and excluding (red solid 
histogram) cluster mergers. The left (right) plot shows the histograms of all Einstein radii larger than 10" (20"). 




-80 -40 40 80 -80 -40 40 80 -80 -40 40 80 



x [arcsec] x [arcsec] x [arcsec] 

Fig. 18: Tangential critical curves of the systems with the largest (left diagram), second largest (middle diagram), and third largest 
(right diagram) Einstein radius in our simulation including cluster mergers. The black cross in the centre indicates the position of 
the main halo. The red crosses mark the positions of all subhaloes that are enclosed in the tangential critical curve. 



Einstein radii to measure the size of tangential critical curves. 
However, modelling cluster mergers by means of simplified toy 
models, we demonstrated that the strict correlation between me- 
dian Einstein radii and lensing cross sections breaks down dur- 
ing these events. In comparison, effective Einstein radii cap- 
ture the evolution of the lensing cross sections notably better. 
Furthermore, we found that the relative orientation of two merg- 
ing triaxial haloes plays a significant role for the evolution of 
the strong-lensing properties. While we were able to reproduce 
an enhancement of the strong-lensing efficiency by factors of a 
few choosing the most beneficial relative orientation, we showed 
that this effect completely vanishes when averaging other align- 
ments. 

The correlation between effective Einstein radii and lensing 
cross sections is tight enough to infer sufficiently accurate es- 
timates of optical depths for giant arcs from cosmological dis- 
tributions of Einstein radii. This approach is computationally 
far less demanding than explicitly calculating individual lens- 
ing cross sections, since the computation of Einstein radii can 
be implemented in a particularly efficient way. We therefore de- 



veloped a new, semi-analytic algorithm for computing cosmo- 
logical distributions of Einstein radii that incorporates the im- 
pact of cluster mergers. Using this new technique, we were able 
to compare the statistical strong-lensing properties of one sam- 
ple of single (isolated) dark matter haloes and a second sample 
of dark matter haloes that includes cluster mergers. Assuming a 
constant source redshift of z s = 2.0, we found that cluster merg- 
ers increase the theoretically expected number of Einstein radii 
above 10" (20") by roughly 36% (74%), indicating that these 
events provide a highly efficient mechanism to enhance the lens- 
ing efficiency of particularly strong gravitational lenses. These 
results clearly show that semi-analytic studies need to take clus- 
ter mergers into account if they aim to question the standard 
cosmological model by comparing the largest observed Einstein 
radii to theoretical expectations. Furthermore, we estimated that 
the optical depth for giant gravitational arcs of those haloes with 
Einstein radii above 10" (20") increases by roughly 45% (85%), 
which highlights the importance of cluster mergers for the statis- 
tics of giant gravitational arcs. Still, we were unable to reproduce 
the doubling of the optical depth found by Fed eli et aL] ( |2006| ). 
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However, this discrepancy can be attributed to the following im- 
portant differences between both studies: 

- Instead of adopting triaxial density profiles with varying 
properties (i.e. concentration and axis ratios), |Fedeli et"aL 



2006 ) modelled clusters as spherically symmetric density 



profiles and elliptically distorted their lensing potentials by 
a constant amount. Moreover, these authors used an alter- 
native concentration-mass relation. Consequently, their lens 
model produces surface mass densities with different prop- 
erties, which are decisive for the strength of the boost of the 
strong-lensing efficiency during cluster mergers. 
Analyses of cluster mergers in numerical simulations reveal 
that (1) the major axes of inf ailing substructures are intrinsi- 
cally aligned with the major axis of the main halo due to its 
tidal field and (2) subhaloes preferentially approach the main 
halo along its maj or axis ( |Lee et al.|2005[|Altay et al.|2006} 
Zhan g et al.|20 09). Motivated by these findings, Fedeli et al 



(2006J throughout modelled mergers by perfectly aligning 
the major axes of both haloes and assuming that the sub- 
halo always approaches the main halo exactly along its ma- 
jor axis. As shown in Section |5^2{ these orientations are most 
favourable and produce the strongest possible enhancement 
of the strong-lensing efficiency during cluster mergers. On 
the other hand, this assumption clearly is an idealization that 
might lead to an overestimation of the net effect of merg- 
ers. Since there is no reliable theoretical framework that al- 
lows us to predict those alignments realistically, we chose 
the fairly conservative approach to sample the relative orien- 
tations of merging haloes randomly, without taking their axis 
correlations and preferred directions of motion into account. 
Hence, we expect that our results slightly underestimate the 
net effect of mergers. 

While we only computed the optical depth o f particularly 
strong lenses with Einstein radii above 10", |Fedeli et~aL 



( 2006 ) determined the optical depth of all clusters at red- 
shifts z > 0.5. Hence, their computation additionally incor- 
porates the following important effect. Cluster mergers tend 
not only to increase the lensing cross section of the strongest 
gravitational lenses, but are also capable of transforming 
clusters to strong gravitational lenses that would otherwise 
be subcritical. More precisely, it may well be that two rela- 
tively small clusters with vanishing individual lensing cross 
sections produce a nonvanishing, joint lensing cross section 
while merging. Consequently, mergers are an effective mech- 
anism to substantially increase the optical depth of clusters 
in the lower mass range. Given the steepness of the mass 
function, these objects dominate the total optical depth and 
thus their contribution is likely significant for the doubling 
of the optical depth found by |Fedeli et al. ( |2006 ). 



Given these remarks, we conclude that our findings do not dis- 
agree with the previous results of Fedeli et al. ( 2006]). 

In this work, we mainly focused on the development of new 
techniques for incorporating the impact of cluster mergers into 
semi-analytic, statistical strong-lensing studies. We did not con- 
duct a detailed comparison between theory and observations, 
in particular with regard to the problem of too large Einstein 
radii. Such analyses will be performed in our follow-up work, 
where we plan to apply methods of extreme value statistics to 
investigate whether or not cluster mergers are required (and suf- 
ficient) to explain the largest observed Einstein radii within the 
fram ework of the standard cosmological model ( Waizmann et al.| 
|20T2] ). 
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Appendix A: Gravitational lensing by triaxial dark matter haloes with generalized coordinate 
transformation 



As mentioned in Section 23 we need to introduce a more general rotation matrix than Qguri et al. ( 2003| ) to describe the coordinate 
transformation between the principal axis frame of triaxial haloes and the rest frame of the observer. Consequently, some algebraic 
expressions in the derivation of the lensing properties of triaxial haloes are more complicated. Adopting the same notation as Oguri 
et al. ( 2003| ), the corresponding expressions are given by 
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Appendix B: A method for accelerating the computation of deflection angles 



Oguri et aL] ([2003) show that the surface mass d ensity profiles of triaxial haloes are ellipsoidal. Let q denote the corresponding 
minor- to-major axis ratio (cf. |Oguri et al.| 2003] Eq. (35)). Then, the components of the deflection an gle a = (a x ,a y ) can be 
expressed as one-dimensional integrals of the convergence /c(£) (see also Schramm| 1990a Keeton 2Q0Ta|: 



a x (x,y) = qxJi(x,y) , 
a y (x,y) = qyJ (x,y) , 



where the integrals J n (x,y) are defined as 



Jn(x,y)= f 
Jo 
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The numerical evaluatio n of these integra ls is computationally quite expensive. Since lim„^o £ (v) = and lim^o k{( ) = oo (cf. 
Eqs. (40) and (41) from Ogu ri et al.| (|2003 )), the integrand is singular at the lower integration boundary. Accordingly the numerical 
integration routines converge only slowly. We found that this inconvenient behaviour can be avoided by two subsequent substitutions. 
First, we invert Eq. (|B.4|> and solve for v: 
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For y t 0, lim v ^o£( v ) - implies that the negative branch must be the correct solution for f «: 1. Since 
f 2 (l - q 2 ) + x 2 +y 2 \ - 4£ 2 y 2 (l - q 2 )\ never vanishes for y + and the solution must be continuous, the negative branch 
is the unique solution for the complete interval. Given this relation, we substitute v by f (v) and obtain 
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The obvious singularity at the lower integration boundary vanishes after this transformation, but the integrand is not smooth at the 
origin yet. We overcome this problem by making the additional substitution f = t 2 , so that 



Jn(x,y)= I 
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These simple integral transformations speed up our computations by a factor of ~ 10. 



Appendix C: Sampling cosmological populations of dark matter haloes 

In several sections of this paper, we needed to sample realistic catalogues of dark matter haloes, where 'realistic' means that they 
should have similar mass spectra as halo populations in numerical simulations and the observed universe. For that purpose, we 
adopted the following c ommon Monte Carlo ap proach (see Oguri & Blandford 2009} for instance), which is based on the Press- 
Schechter mass function (Press & Schechter|1974|). Let us assume we need to draw a random population of dark matter haloes within 
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the mass range [M m i n , M max ] and the redshift interval [z m in, £maxL Then, we first subdivide both ranges into smaller, equidistant bins, 
adopting typical bin sizes of Az = 0.02 and A(log(M)) = 0.02. Next, we calculate the mean expected number N of haloes in each 
bin using the Press-Schechter mass function dn(M, z)/dM, 

d 2 N dn(M,z) dV 

dz dM dM dz 

where dV/dz denotes the differential comoving volume. After that, we generate a random integer number N from the Poisson 
distribution with mean N. Finally, we sample N haloes within the current bin and assign them uniformly distributed masses and 
redshifts out of the corresponding ranges. 

This Monte Carlo method is simple and efficient. It reliably generates random halo catalogues whose deviations from the Press- 
Schechter mass function remain within the limits of Possion noise. 
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